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Abstract: We consider the application of the polyharmonic subdivision wavelets 
to Image Processing, in particular to Astronomical Images. The results show an 
essential advantage over some standard multivariate wavelets and a potential for 
better compression. 
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1 Introduction 

In [5] we have constructed the basic elements of the polyharmonic Wavelet Analysis. 
There we have provided all details of the construction of the filters for the whole 
infinite family of father and mother wavelets. Here we will finalize the multivariate 
construction and provide experimental evidence of the power of the new wavelets. 

2 Polyharmonic Wavelet Analysis 

The Polyharmonic Wavelet Analysis arises in the context of the polyharmonic sub- 
division. In general, it is obtained when the (one-dimensional) polynomials of degree 
2N — 1 are replaced by polyharmonic functions of degree N. The main element of 
the construction is the application of the Fourier transform in n variables in M""*"^. 
For simplicity we assume that the Image which we consider is a function u (t, y) , for 
i G M and 27r— periodic in e M". We have 

u{t.y)= {t) e^^^'^'^ where the Fourier coefficients are (1) 

r,eZ" 

For every fixed rj G Z"^ and ^ = \ri\ we make Wavelet Analysis which is based on 
the non-stationary father and mother wavelets {v'rn (^)}„i>o ^^^"^ {'^mit)}„i£Z '^^ 
Theorem 4. Recall that they have refinement coefficients Afj™'"'' and (—1)"' Mj™|^'^ 

respectively, where the polynomial M^™!''' (2) was defined at the end of 5 . This 
gives the wavelet expansion 

u (t.y) = ^v^m^A'J it - j) e'^""^"^ ■ (2) 



We summarize the results in the following theorem, for the proofs we refer to [3] . 



Theorem 1 For every ^ > the functions V'm {t) generate a non- stationary MRA 
of L2 (M) by means of the definitions 

Vi := dosL,(R) {ifi, {t - i) -.jeZ} 




L2 (M) = doSi2(R) I 

Here cZos^^CR) denotes the usual closure in a linear and topological sense with respect 
to the space L2 (M) . Respectively, the Wavelet Analysis spaces Wm defined by means 
of 

are generated as 

^ closL.iM) {^i it-j):jeZ}. 
As in the polyspline Wavelet Analysis studied in 4 , the one-dimensional MRA 
y'''' which is generated for every 77 S Z" is related to a " polyharmonic MRA" PVj 

as PV, = e,,,„ yl"! by means of the Fourier transform formula |. Respectively, 
we have also the polyharmonic Wavelet Analysis defined by PVj ^ PWj — PVj^i 
and reduced to infinitely many one-dimensional Wavelet Analyses for every 77 g Z" 
in the form PWj = ^j''' means of the Fourier transform formula 

Respectively, the space 

L2,per ^{f{t,y):teR, yE M"; / is 2^ - periodic in y} 

is generated by an infinite amount of orthogonal (non-stationary) mother wavelets 

(Oe'^"'^^ :r/GZ", tti e z} . (3) 

Although these wavelets do not have a compact support in direction y, it is important 
to note that they have "elongated support" in the y direction, which plays an essential 
role for the sparseness of the "edge representation" in the images, a point of view 
nicely described in [2], see also [B]. 

3 The algorithm 

From the above it follows that the algorithm for Decomposition and Reconstruction 
of the Polyharmonic Subdivision Wavelet Analysis (PhSdWav) consists of 
a Fourier transform in direction y and consecutive application, for every rj g Z", of 
one-dimensional Wavelet Analysis in direction t. It consists of the following steps: 

1. As usually we take as first approximation to the image the image itself where 
the expansion is in terms of the scaling functions I ipl^i (t) I 

2. Transforming to Wavelet Domain - this is the main difference with the usual 
approach. Roughly speaking, we obtain (t) in formula ([T]) and make its 
Wavelet Analysis for every rj G which gives us the coefficients ^,-i,m.j in the 
wavelet expansion ([2]) . Let us note that these coefficients are in general complex 
numbers. 

3. Filtering the small (insignificant) coefficients by using a threshold parameter. 

4. Quantizing the remaining coefficients. 

5. Encoding by arithmetic coder. 
The algorithm was written in Matlab. 
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4 Experiments with images 



We apply our Wavelet Analysis to different classes of images - synthetic and real. 
We compare the results of our method PhSdWav with parameter N = 9 (see [S]) 
against the 2c? wavelet transform with Daubechies wavelets db9 (which are tensor 
product wavelets and correspond to the one-dimensional Daubechies wavelets with 
parameter iV = 9). We use the standard function dwt2 in Matlab which contains a 
large class of 2d wavelets, and db9 is one of the most used representatives. In all 
Figures below we first display on top-left the original image, on top-right is the 
result of PhSdWav compressed representation obtained by wavelet thresholding, and 
on the bottom is the result of the db9 compressed representation obtained by wavelet 
thresholding. In order to make comparison between our method PhSdWav and db9 
we perform thresholding which results in equal PSNR. 

4.1 Elongated support tests on 64 x 64 pixel edges 

For the proper compression and representation of images it is of fundamental impor- 
tance to understand how the method applies to " edges" . For that reason we will first 
show the performance of our method to some synthetic edge images. 

4.1.1 Vertical edge 

In this test the image is a vertical edge which is a function f^ert {t, y) defined in 
the rectangle [0, 1] x [0, 1] by 

Jo for < t < 1/2 
hert{t,y) forl/2<t<l- 

We take a sampling of / {t, y) by a 64 x 64 pixel image. The performance of the 
polyharmonic Wavelets and the db9 wavelets is provided in Table 1 and Figure 1 
below. 



Table 1. 





Num. cocffs 


PSNR 


Entropy 


Pol. Wav. 


48 


120.7379 





db9 


960 


121.7412 


l.Ole+003 



As we said above the wavelets of our basis (|3| show elongation in direction y and 
this makes the above image fvert perfect for representation and compression with 
the polyharmonic subdivision wavelets. We see that the quality of representation by 
means of PhSdWav (measured by the PSNR) is the same as the quality of repre- 
sentation by means of the standard 2— dimensional Daubechies wavelet db9 but the 
number of coefficients is 20 times less! 

The image fvert would correspond to the one-dimensional image provided by the 
Heaviside function, / (i) = for < t < 1/2 , / (i) = 1 for 1/2 < < < 1. The wavelets 
have an excellent performance for such images. 

Just the opposite is the next example. 

4.1.2 Horizontal edge 

In this test the image is a horizontal edge, which is a function fhor (t, y) defined in 
the rectangle [0, 1] x [0, 1] by 



fhor {t, y) := 




for < 2/ < 1/2 
for 1/2 < y < 1 ■ 
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Figure 1: Vertical Edge 



We take a sampling of fhor {t, y) by a 64 x 64 pixel image. The performance of the 
polyharmonic Wavelets and the db9 wavelets is provided in Table 2 and Figure 2 
below. 



Table 2. 





Num. coeffs 


PSNR 


Entropy 


Pol. Wav. 


2560 


121.4740 


2.528+003 


db9 


960 


121.7412 


l.Ole+003 



We see now how inefficient the PhSdWav are for the representation of fhor com- 
pared to db9 wavelets, since for the same PSNR we need almost 3 times more coeffi- 
cients. This result is expected since PhSdWav has a very strong vertical orientation. 

4.1.3 Skewed edge 

Somewhere in the middle is the case of a skewed edge. In this test the image of 
skewed edge, which is a function fskew {t^u) defined in the rectangle [0, 1] x [0, 1] 
by (this is an approximate definition) 

. ity)-[^ for0<i<i(y + 2.5) 

We take a sampling of fskew {t, y) by a 64 x 64 pixel image. The performance of the 
polyharmonic Wavelets and the db9 wavelets is provided in Table 3 and Figure 3 
below. 
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Figure 2: Horizontal Edge 



Tabic 3. 





Num. cocfFs 


PSNR 


Entropy 


Pol. Wav. 


512 


126.8204 


386.7158 


db9 


940 


121.9320 


977.3171 



We see that although the edge is prettily skew (the slope is approximately tan a = 
6 ) the PhSdWav have a very nice performance compared to db9. 

4.2 Lena image 

Let us make a more detailed analysis of the seminal image of Lena in 128 x 128 pixels. 
We provide the Compression Ratio and the Bits Encoded for both methods in Table 
4 and Figure 4 below. 



Table 4. 





N. coeffs 


PSNR 


Entropy 


Cmpr. rat. 


Bits Enc. 


Pol. Wav. 


4,150 


77.7338 


5.15e+004 


5.9815 


42,802 


db9 


5,187 


77.6742 


5.65e+004 


5.9562 


44,012 



We see that for the same PSNR our method PhSdWav is better - it needs 20% 
less coefficients. This is an astonishingly good result since as we have seen in the 
experiments with the synthetic edges in Section |4.1[ our method of PhSdWav is 
adapted to the vertical direction y and the Lena image does not have a dominating 
orientation of the edges. 

4.3 Astronomical Images 

Hereafter we analyze different types of astronomical images. 
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Figure 3: Skewed Edge 

4.3.1 The Pleiades star cluster, a part with faint objects only, sized 128 x 
128 

First we consider an image which is a part of the Pleiades image containing only faint 
stars. The results are provided in Table 5 and Figure 5 below. 



Tal)l<- 5. 





N. coeffs 


PSNR 


Entropy 


Cmpr. rat. 


Bits enc. 


Pol. Wav. 


2,663 


77.5805 


2.93e+004 


7.6028 


33,456 


db9 


5,187 


76.3970 


4.52e+004 


6.9317 


37,818 



Wc sec that for (almost) the same PSNR wc have much less coefficients, the 
compression ratio is higher, and our entropy is much less, which shows a better 
potential for compression. 

4.3.2 The Pleiades star cluster, a part with bright objects, sized 256 x 256 

Next we consider a part of the Pleiades containing several bright objects. The results 
are provided in Table 6 and Figure 6 below 



Table 6. 





N. coeffs 


PSNR 


Entropy 


Cmpr. rat. 


Bits cnc. 


Pol. Wav. 


15,159 


50.0717 


2.89e+005 


13.8338 


73,750 


db9 


18,502 


49.4214 


3.57e+005 


8.5379 


122,814 



Again wc sec that for the same PSNR we have not only less coefBcients but also 
a much higher compression ratio. 
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Figure 4: Lena Image 
4.3.3 The Pleiades star cluster, chain, sized 128 x 128 

Here we consider a chain image of a part of the Pleiades image. The results are 
provided in Table 7 and Figure 7 below 



Table 7. 





N. coeffs 


PSNR 


Entropy 


Cmpr. rat. 


Bits enc. 


Pol. Wav. 


2,764 


88.6818 


2.48496+004 


8.8706 


28,528 


db9 


5,340 


87.5730 


5.7468e+004 


8.4636 


30,973 



For the same PSNR, the number of coefficients of our method PhSdWav is twice 
less than the number of db9 coefficients. This result is expected since the image is 
vertically oriented. However the compression ratio is not much bigger apparently 
due to more noise in the background of the image. 
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Figure 5: Pleiades star cluster, faint objects, 128x128 pixels 
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Figure 6: Pleiades star cluster, bright objects, 256x256 pixels 
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Figure 7: Pleiades star cluster, chain image, 128x128 pixels 
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